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Abstract 


A computer program aimed at the phase separation between gas and liquid at 
zero gravity, induced by vortex motion, is developed. It utilizes an explicit 
solution method for a set of equations describing rotating gas-liquid flows. 

The vortex motion is established by a tangential fluid injection. A Lax- 
Wendroff two-step (McCormack's) numerical scheme is used. This program can be 
used to study the fluid dynamical behavior of the rotational two-phase fluids in 
a cylindrical tank. It provides a quick/easy sensitivity test on various 
parameters and thus provides the guidance for the design and use of actual 
physical systems for handling two-phase fluids. 

Key Words: computer code; gas-liquid separation; numerical modeling; two-phase 

vortex motions 
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I. Introduction 


Mechanical systems have been devised for producing artificial gravit 
fields to spin-up liquids in containers. These involve rotating mechani 
are cumbersome and, more importantly, have moving parts that wear out. 

Here, liquid rotation created by fluid injection is consider ?d, 'he detail 
analysis of the two-phase vortex model can be found elsewhere [1]. In this 
report, the details of the computer code are described. 

The computer program was developed to study the fluid dynamical behavior of 
two-phase fluids in a tank at zero gravity. The phase separation between gas 
and liquid, induced by vortex motions, is of primary interest. The program 
utilizes an explicit solution method for a set of equations describing rotating 
gas-liquid flows. The vortex motion is established by a tangential fluid 
injection. A Lax-Wendroff two-step (McCormack’s) numerical *scheme is used in 
the computer program. This scheme uses a conservation form of a system of 
equations together with an auto time step feature. 

The program was developed and tested on an HP-1000 minicomputer. The HP- 
1000 ’s FORTRAN 77 is based on the American National Standards Institute (ANSI) 

77 standard programming language FORTRAN (ANSI X3. 9-1978). The HP FORTRAN 77 
has extensions to provide a more structured approacl” to program development and 
more flexibility in computing for scientific applications. It fully implements 
the Military Standard Definition (MIL-STD-1753) of extensions to the ANSI 77 
standard. In order to make the computer code more useful for other computer 
systems, modifications have been made so that the code is closer to the ANSI 77 
standard and thus less system dependent. Some limited extensions are still kept 
in order to produce the code in the HP-1000. Since the graphic routines are 
system dependent and must be modified with their equivalents at each computing 
facility, the original graphic code has not been included in this report. All 



lines preceded by "“V" are originally adopted to use the vector operation 
package supplied by Hewlett Packard. The speed of the code can be increased by 
replac? ig many "do-loop” operations in the code with high speed vector 
operations. Little effort is required to incorporate the vector operation into 
the code if the vector operation package is in the system. 

Thus with limited effort, the program can be adapted easily to most systems 
accepting the ANSI 77 standard FORTRAN. For example, the EMA (Extended Memory 
Area) statements may have to be removed from the code fcr some computers. Also 
double precision real numbers (Real * 8) could be replaced by single precision 
real numbers. 


II. Model Equations 

The vortex induced model is based on a two-phase, two-fluid continuum [2]. 
It incorporates several interactions between phases; namely fluid drag and 
virtual mass effects and it can be modified to include additional interaction 
effects. Detailed analysis of the model has been reported in Ref. 1. A brief 
summary of the system of equations Is given below. 

The equations for the conservation of mass and momentum for the two fluid 
two-phase model in an one-dimensional, axisymmetric case 

(i * e * h • Is • 0) are: 


i 


3ro k 

at 


3r Vrk 

3r 


0 


2 



9ra. V Bra. V „ „ 

k rk . k rk „ 2 „ _ 3p 

— — + — - a V » -a, C . r -r*- 

3t 3r k 9 k k pk 3r 


2 3r VrrZ 

+ a k ^ C kZ ( dr VeeZ + Vz B rZ r) 


“k C dk r (V rl " V r2> 


3ra k 0k 3ra k V rk V 0k v u r r / 3ra Z T r0Z - r \ 

~lt + 3? + a k V rk V 0k ' a k EC ki ( ~ a Z T r0Z + Vz B 0Z r) 

Z»1 


for k » 1 and 2 and with 


+ a k c dk r (V 91 - V 


c pl * (a lV2 + A a )/<pt> 


C p2 “ (a lVl + A a )/<p > 


C 11 * (a 2 P 2 + V /<p > 


C 12 ■ C 21 " A a /<P > 


C 22 ' ( Vl * A a )/<P > 


n 



~dl 


= -a. 


./<P 2 > 


'62 


ViV <p > 


and 


<P 2 > 


a iVi p ? * 4 * 


'Vi 


« 2 p 2 ) 


The effective stresses are modeled as 


t , » 2p, 3V , /3r 
rrk ^k rk 


’rOk ' T erk ■ “k r 3 < V 9 k / '-> /3 '- 


r ani « 2 y, V , /r 
00k k rk 


with 



and the interfacial forces are modeled in the form of 


VWV * 4 a at - V- 

H. is the force density acting on the phase 1 by the phase 2. A 
added mass and drag coefficients, respectively. 


The incompressibility condition is reduced to + a 2 '/ p2 


and A . are the 


■ Q /r, where 
r 


Q is tne net radial outflow, 
r 

In the program Q p - 0 is assumed, since the mixture pumped out is injected 
immediately back into the tank at the nearby location. The net volume or mass 
in the system is effectively unchanged e A cept for the net change on the angular 
momentum. Thus, the pump system (withdrawal and injection) acts as a body force 
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on the mixture at the nozzle location. The net momentum gain is thus the 
momentum introduced into the system minus the local momentum pumped out. 
Therefore, we will model this pumping dynamic by body forces without considering 
the mass transfer. That is, the body force density a^p^B^r will be replaced by 

Vk v i - - 

the net momentum gain, — r — - (V.n - l r ) at the nozzle location, where V is the 

J J 

injection speed. 

III. Numerical Method 

The complete solution of the complicated system of equations can only be 

obtained through numerical methods. An improved Lax-Wendrof f , two-step scheme, 

(also referred to as MacCormack’s method) [3, 4] is adopted for solving this 

time-dependent problem. This non-centerecf differencing scheme, using a full 

step bacKward prediction and forward correction version, requires no explicit 

artificial viscosity if a proper stability condition is satisfied, Using this 

technique for solving fluid flow problems is very efficient and has been in 

widespread and successful use for some time. It is good both for the time- 

accurate computation of steady and unsteady flow problems. The general features 

of the scheme are: i) its explicitly conservative form, ii) it is a two-step 

predictor-correction type, iii) it is three point, two level - that is, the 

solution of f^*^ - at level n+1 depends only on three values of f^ 1 at level n, and 

lv) it is second-order accurate in time and in space. 

For using the MacCormack’s numerical technique, the system of equations can 

be expressed in the conservative form as: 

W - F + P ♦ gG ♦ S 
t r r r 

Here the subscripts (t and r) denote partial differentiation witn respect to t 
and r, respectively, and W, F, P , gG^ aid S are column matrices with five 
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i 


/ 


elements. Ail the components of F, P , «G ana S can be regarded as functions 

r r 

of the components of W which are the independent variables. The fundamental 
theory of the MacCormack's scheme is briefly given below. 

For second order accuracy, the solution could be written as 

W 1 - W° ♦ At W 0 * W 0 

t C. L u 



At W_°) 

U v 


- ± ( W° 


At W t °) 


1 

2 


(W^ 


At w p ) 


where 


and 


- | (W p ♦ W c ) 

W p ■ W° * At W 0 is the predicted value, 

u 

W c - W° ♦ At W t p is the corrected value. 


The superscripts denote the time-level of the information and subscripts denote 
the partial derivative with respect to either time t or space r. Specifically, 
superscripts 0 and 1 are the initial and the completely advanced time (here two 
steps) plane; p and c are the predicted (1st step) and corrected (2nd step) time 
plane. Thus, W^ 0 is the time derivative of W evaluated at the initial time, and 
W t P is time derivative of W evaluated at the predicted time. 

Fig. 1 shows the diagram of the two step difference scheme used in the 
computer program. Due to the difference scheme, the spatial location after each 
step in time is a half grid off from the original one. Thus, the spatial offset 
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which resulted from a backward predicting step will cancel with those of the 
forward correcting step. 

Numerically, the predicted values are 


wP „ I ( W P + u p ) 

i 2 ' "i-1/2 "i+1/2' 


where 


*1-1/2 ■ I <W 1-1 * W l°> * “ "t 


1 {h 0 + m °) ♦ it r ( F 0 - F ° ) + — iili- (G 0 - G 0 )] 

2 1 i-1 i ' Ar i i-r 2 v i i-l' J 


+ AL ( sP + S ° ) + At P P 

2 ts i i-l } At i-1/2 . 


and the corrected value is evaluated at the predicted time place, that is at 


l + l/2‘ 


Thus 


wr - 


W 0 ♦ At W P 

1 v 


. u 0 + At r(F P 
W i Ar L(F i+l/2 


( g P ♦ g P ) 

P ^ g l+l/2 S l-l/2 ; p p }] 

F i-l/2 ) 2 ^ 1+1/2 G i-1/2 ;J 


+ 


At 

2 


(S 


i+1/2 


+ S v ) 
S i-l/2^ 


+ At P 


c 

i 


Here P p /2 and P. c are the pressure correction terms at half and full time 
steps respectively. Thus, for each time step, the advance is carried out in two 
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steps: a full step backward predictor, and then a forward corrector. As 

indicated in the diagram, the subscript i is the regular mesh spatial location 
at which solution is to be advanced, i ♦ 1 is the spatial location of regular 
mesh points immediately to the right and left of the location i, i ♦ 1/2 is the 
location midway between I and i ♦ 1 or between i - 1 and i at the predictor 
plane. Thus, for each time step as the procedure advanced, the outermost data 
points at the boundary are not updated through the numerical scheme. The values 
at the boundary are to be given through some suitable boundary conditions. The 
numerical procedure utilizes a uniformly preselected spatial mesh and variable 
time increment. To avoid a singularity at the center of the core region, a 

finite radius is used for the inner boundary. The tank radius R Is the outer 

boundary. The time step is determined at each time step to ensure numerical 

stability [5]. For a finite grid size Ar, the maximum* time step At is given by 

At k “ 1/£ -I C dJ + l V rJ /Ar * ~~2 (o l u l C kl * ®2 W 2 C k2^ 

Ar 

where k - 1 and 2. The minlmimi At^ (with some rounding off) is used for the 
time step. Normally, the technique with the time step condition gives fairly 
good numerical stability. However, in critical conditions nunerical damping can 
be added either for damping oscillations due to large gradients or for 
accelerating the calculation by increasing the time step. A damping factor, 0 
thus was added in the program as 

«l i0 - Hj 1 (1-D) ♦ - wj) D 
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X ID 

where is the value obtained based on the two-step scheme, and W A is the 

value after the damping factor D Is added. A typical value of D - 0.2 can be 
used for debugging the program. If no damping factor is desired, D - 0 should 
be used. 

The computer program was written in a Fortran 77 based computer code. The 
code will permit evaluation of the effects of various parameters which control 
the fluid dynamical behavior. These include tank size, fluid properties, such 
as density and viscosity, etc., characteristic gas bubble and liquid drops 
sizes, and relative location of injection nozzles. 

A sample input and its output are shown on Exhibits A and B, respectively. 
The initial conditions for the gas and liquid volume fractions are taken to be 
25$ gas and 75$ liquid. These fractions are uniformly distributed over the 
circular cross-section of the cylindrical tan*. Initially both fluids are at 
rest. Other parameters can be fouhd in the sample input In Exhibit A. The 
resulting velocity distributions and gas volume fraction as function of time for 
the sample run are given in Figs. (2) and (3) respectively. The velocity 
distributions are displayed along equally spaced rays at different times to 
enable clear observation. These velocity vector fields indicate all flows are 
primarily in angular rotation with gas phase tending to move inward and liquid 
phase trying to move outward, as expected. As the result of these radial 
movements, the volume fraction distribution is also changed with time. And as 
expected, the gas volume fraction is Increasing at the inner region and 
decreasing at the outer region as shown in Figure 3* More detailed results have 
been reported in Ref. 1 . 





IV. Program Details 


The complete computer code is listed in Appendix A. The code consists of a 
main program, GLVM and several subroutines. It is written in subroutine form 
such that each subroutine performs an individual task. Each logical part is 
clearly Isolated and it can be easily modified to reflect different modelings 
for the interfacial forces. The interactive input mode with self-instruction is 
used for easy parameter insertion. Many instructive internal documentations are 
included in the program. In the following, each subroutine is listed with a 
brief description of its major function. 

1) GLVM, the Main Program. 

•To initialize data and start the program: Logical unit to save data (LUS), 
(Logical Unit is 1 for terminal, and 6 for printer), job identification 
notes (NOTES), data file name for saving data (NAMR), initial time (TO), 
final time (TMAX), time interval for data output (DTPRT), etc. 

•To oontrol the calling sequences to the other subroutines. 

•To check the time step. 

•To save, print (and plot) the output data. 

•To obtain the predicted and corrected values in the two step, numerical 
scheme. 

•To impose boundary condition. 

•To update the data, time, and step number for the new time step. 

•To provide a shutdown procedure either in normal (e.g., t > t ) or 

UlftA 

abnormal (e.g., At is too small) conditions. 

?.) INIT, Initialization. 

•To input the test parameters, initial conditions and set-up the initial 
column matrix N. 
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Default values are provided for most of the parameters. The default values 
are listed at eacii interactive input step. If the default value is acceptable, 


a comma " 
Some 

ALMT 

DAMP 

DEN1D2 

DO 

DS 

EVF 

GAMMA 

IVTX 

« 

IW 


MM 

.J 

MLEF 

MU1D2 

NA, ND 


" is inputted. 

of the relevant symbols used are listed below: 

Limit values of a^, ALMT(l) < < ALMT(2). 

Numerical damping factor. Normally set to 0. 

Density ratio, P-^P 2 

Y 

Base diameter, i.e., d k ■ d Qk a k 
Density scale - p 2 

u^/u effective eddy viscosity factor. 

Diameter exponent Y 

Type of simple initial flow: 0 - at rest, 1 - simple rotation, 2 » 
Hammel-Oseen Vortex, 3 - G.I. Taylor Vortex. This is needed only when 
there is no data file (NAMR) given for an initial condition. 

Boundary condi. ion at wall (for tangential component). 0 - free-slip , 
(no-skin friction), 1 - non-slip. Also, when 1^ ■ 1, a factor of (1- 

r)®*~ was included on IVTX flow to simulate an Initial power law 
boundary layer. 

Size of data arrays. Mil > NG. MM appears in many subroutines. 
u k , dynamic viscosity for phase k. 

(1 + EVF) * MU. Effective viscosity. 

Pj/Ug. Viscosity ratio. 

n ,n , Exponents for weighting function for drag and added mass 

3 0 

coefficients, and A.. 

a a 
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NAMR Data file name for initial condition, if any. 

NC Number of grid points used. 

OMEGA w, initial rotation speed, V 0 - «r, when IV7X > 0. 

PS Pressure scale, DS * VS^ 

QJ,VJ Injection flow rate and speed. 

RE Reynolds number, p^V^R/^ 

RJ1.RJ2 Jet opening, RJ1 < r < RJ2 

RPEAK Location of the peak speed of the initial vortex, if IVTX > 1 

RTANK Tank radius - Length scale LS. 

VPEAK Peak speed of the initial vortex, if IVTX > 1. 

TJ1.TJ2 Injection time, TJ1 < t < TJ2. 

TS Time scale - L S/VS 

VS Velocity scale. 

. The format of the data file for the initial condition (if any) is a six 
column and NG row data file, where NG is the number of grid points. The column 
sequence is K, ( K ) , V rl (K), V r2 (K), V 01 (X), V 02 (K), where K is the grid point 

number, and the rest of the terms are the gas voliae fraction, radial gas 
velocity, radial liquid velocity, tangential gas velocity and tangential liquid 
velocity respectively. The data format is free. 

3) DELA, AA 

To determine the fraction of grid site in which the injection is made. 

0 < AA < 1. This is used to define the location of Jet. The region of 
injection could cover several full or fractions of grid sizes. 

*) DERIV1, Derivative 

To get the first derivative of a data array using a center difference scheme 
except the two end points in which three points near the boundary are used. 
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5) DGCOEF, Generalized Coefficients 

To calculate the added mass, drag and all the generalized coefficients (A 

d 


A^ and C^). This is the heart of the modeling. 
The effective coefficients are modeled as: 


A a ■ A al“al * A a2 u a2 


A d " A dl“dl * A d2 w d2 


A al * a 1 a 2 p 2 / ^ a l + 2a 2 / ^ 1 + 


A a2 " “ 1 « 2 p i / ^ 01 2 + 2ot i / ^ 1 + 3a 2 )) 


A dl ' 18 u 2°l /d l “2 


A d2 - 18 lyij/Ul - a 2 /0.8: 2 d 2 2 ) 


na / / na na * 

“al ' “2 /( “l * “2 > 


w I m 1 — w 

a2 al 


nd,/ nd ^ nd. 
“dl ’ “2 /( “l 4 °2 ! 


d2 A dl 

6) DWPDE, Partial Differential Equations. 

To evaluate the values of the increments on the column matrix W from the 
partial differential equations. This is the major part of the McCormack's 
scheme. In each complete time step this routine will have to be called twice 


7) FNDDT, at 

To determine the suitable time-step size. 

8) FSOFW, Column matrices F and S. 

To determine the convective matrix F and the source matrix S. 
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9) JET, Injection. 


i 


To determine the momentum source due to the jet injection. 

10) SIZES 

To determine the gas bubble and liquid droplet sizes. In the model the 
sizes were modeled to be functions only of the volume fraction, i.e. 

"k • d ok \ y - 

Different models for size distributions could be easily adopted here. 

11) TAUOFW 

To determine the stress tensor t and its derivative. 

12) UOFW 

To convert the column matrix W into the physical independent variables, 
such as a, V , V a . 

V. Summary 

A computer program aimed at the phase separation between gas and liquid at 
zero gravity, induced by vertex motion, is developed. The vortex motion is 
created by fluid injections. The computer program uses a FORTRAN 77 based code 
and HP-1000 minicomputer. It is flexible and accepts various input parameters 
for different flow conditions. Other Interaction effects can also be added or 
modified easily. This program can be used to study the fluid dynamical behavior 
of the rotational two-phase fluids in a cylindrical tank. It provides a 
quick/easy sensitivity test on various parameters and thus provides the guidance 
for the design and use of actual physical systems for handling two-phase fluids 
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Code Listing 


T-00004 IS ON CR T4 USING 00126 BLKS R-0000 
FTN77 

•EH* /DAT*/ . /UWM/ , /COEFF/ , /SOURCE/ , /FANDS/ , /TAU/ 

•FILES 1 .2 

PROGRAM CL VM(,99>, <860425. 1537) 

C 

C THIS PROGRAM UAS DEVELOPED TO STUDY THE FLUID DYNAMICAL BEHAVIOR 
C OF A ROTATIONAL TWO-PHASE FLUIDS < GAS/L IQUID ) IN A CYLINDRICAL TANK. 
C THE VORTEX MOTIONS ARE ESTABLISHED BY TANGENTIAL FLUID INJECTION. 

C THE PROGRAM WAS DEVELOPED ORIGINALLY BY T.T. YEN Of- NBS 
C IT WAS BASED ON HP'S FORTRAN 77<ANSI 77+MIL-STD-l 733 ) 

C 

C WHEN WHO WHAT 

C 8302XX TTY ZE O-G FUEL TRANSFER, START-UP STAGE. 


C LAa WCNDROFF 2-STEP SCHEME (TULL STEP PREDICTION+CORRECTION ) 

C WITH NUMERICAL DAMPING FACTOR < NORMALLY SET TO ZERO) 

C CONSERVATION FORM, VAR I ABLE < AUTO ) TIME STEP 

C REAL»0 

C INTERFACIAL FORCES: DR AC, ADDED MASS 

C PUMP CONDITION: MOMENTUM SOURCE BUT NO MASS SOURCE 

C 850713 TTY GENERALIZED EQUATIONS AND COEFF . Cij 
C 851018 TTY IN ANSI 77 STANDARD < WITH A LITTLE EXCEPTION FOR 
C TESTING IN HP-1000) 

C 

C ***** INTERNAL SUBROUTINES ***** 

C DELA, DERIV1 , DGCOEF, DWPDE , FNDDT , P SOFU , 

C INIT, JET, SIZES, TAUOFW jnd UOFW 


C *** MOST OF /HE LIST OF NOTATIONS ARE GIVEN IN SUBROUTINE XNIT 
C 

CHARACTER NAMR*t 6 , N0TES*72 
INTEGER I , I OS , J , JT I ME < 3 ) , K , MM , NPRT , NT 
Z , IU, LUP , LUS , NG , NGM1 ,NCM2 

PARAMETER (MM-101) 

R EAL*6 BA\ MM ,2 ) , DDT , DT , PTMAX , DTMIN , DTPRT , PZERO 
X , DW1 , T , TMAX , TPRT, VJT, VDR<2) 

Y ,RJ1 ,R J2 , TJ1 , T J2 , QJ , VJ 

1 , ALMT ,U,V,ALP,P,R, U , UP , WN , DW, RDP ,RH, F,S 

4 , BR , BRH, RHO , MUEF , V 1 8 , NA , ND 

5 , DO, GAMMA, DAMP, DR 

6 , TRR , TRA, TAA,RTRR ,RTRA, C, CPA, CD 

COMMON 

Y /JETS/ RJ1 ,RJ2,TJt ,TJ2,QJ,VJ 

Z /CONTP/ IW, LUP ,LUS , NG , NGN! , DAMP , DR 

1 /ALPLMT / ALMT <2 ) 

2 /COEFF/ C<MM,2,2),CPA<MM,2),Cn<MM,L) 

3 /DATA/ U(MM,2),V(MM,2), ALP < MN , 2 ) , P < MM > , R < MM > 

4 /DRAG1/ R HO < 4 ) , MUEF (2) , V 1 8 < 2 ) , NA , ND 

6 /FANDS/ F(MM,3),S<MM,5> 

7 /90URCE/ BR <MM) , BRH < MM ) 

8 /TAU/ TRR < MM, 2) ,TRA<MM,2) , TAA<MM ,2 > , RTRR <MM, 2 > , R7RA<MM , 2 > 

9 /WUWW/ W < MM , 3 ) , WP < MM , 3 ) , WW < MM , 3 ) , D W < MM , 3 ) , R DP < MM ) , R H C MM > 

EQUIVALENCE (UN, BA) 


C »•*••• RHOU) < RHO< 2) <i.e. PHASE-1 “GAS , * > HASE-2»LIOUID ) 
C DTPRT TIME STEP FOR PR INTOUT < AND PLOT) 

C 
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003? 

0060 

0061 

0062 

0063 

0064 
0 065 
0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 
0079 

0079 

0080 
0091 
0082 
0083 
5004 
0083 
0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0 096 

0097 

0098 
5099 
0100 
0101 
0102 

0103 

0104 

0105 

0106 

0107 

0108 

0109 

0110 
0111 
0112 

0113 

0114 
0113 
0116 

01 17 
0118 


LUP-1 I LU FOR PRINTING DEBUG DATA<*1 TERMINAL) 

LUS-6 » LU FOR STORING DATA< *6 PRINTER) 

7 F0RMAT(2X,A > 3I5) 

8 F0RMAT(2X,A,3( 1PE12, 4) ) 

WRITE<1,7> 'Enter lu far saving data. D.F. *',LUS 
READ(1,*> LUS 

C KEET JOB TIME FOR FUTURE REFFERENCE 

C CALL EXEC <11, JTIME , JTIME ( 1 ) ) 

IF (LUS , NE . 1 .AND. LUS . NE . 6) 1 HEN 
C «a****** Define a file name for strcing output a#*****##*#***## 
URITE<1,'(2A) ’> 'Enter FILE NAME for saving data.' 

READ< 1 , ' ( A) ' ) NAH« 

LUS-90 

OPEN < LUS ,FILE*NAMR , IOSTA1 -IDS , STATUS* ' NEW ' , ERR*9°9 > 

ENDIF 

WR I TE ( 1 , ' < A ) ' > ' Enter N0TES(<73 CHAR.) for The job' 

READ( 1 , ' < A) ' ) NOTES 
WR ITE < LUS , ' < 3H1 ,A>') NOTES 

C WR I TE ( LUS , '(514)') JTIME 

C To set-up the initial condition. 

CALL INIT 
'NGM2-NG-2 

VDR < 1 ) *2 . »MUEF ( 1 ) /DR*#2 1 For determining time step 

VDR(2)-2. »MUEF<2)/DR«*2 

T-0 . ! INITIAL TIME 

TMAX-3. 

DTPRT-0.2 

UR ITE< 1,8? 'Enter INITIAL and FINAL TIMEs. D.F.*' , T,TMAX 
R£AD<1,*> T , TMAX 

WRITE ( 1 , 8 ) 'Enter TIME STEP for output. D.F.*', DTPRT 
READ( 1 , ») DTPRT 

DTMIN-1.0D-6 » SET MINIMUM TIME STEP 

DTMAX-DTPRT 

NT*0 \ Time step number 

TPRT'T 

NPRT*0 

PZERO*0.D0 f Pressure at centtr core 

DT*DTMIN 

10 NT*NT*1 

CALL UOFW(U,R,*JG> 

C PRINT OUT AT SELECTED TIME 

IF ( T . GE . TPR X -OR. T . GT . thpx) THEN 
IF (NT .GT. 1) THEN 
NPRT-NPRTM 

TPRT-TPRT+DTPRT»< 1 , *DNJNT < < T-TPR T ) /DTPRT > ) 

21 FORMAT < 1 HI , 2< A5 , 14 , AS , IFF .3)) 

WR ITE ( LUP ,21 ) 'NP-' , NPRT , ' T* ' , T , ' NT- ' , NT „ ' DT* ' , DT 
W*ITE<LU3,21 ) 'NP*' ,NPRT, ' T- ' , T , ' NT* ' , NT , ' DT- ' , DT 

WRITE (LUS, ' <A3,A6,5A1 1 ) ' > ' J ' , ' ALP 1 ' , ' Ut ' , ' U2' , ' U1 ' , ' U2 ' , ' P ' 
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f 


f 

i 



01 19 
0120 
0121 
0122 

0123 

0124 
01 25 
0126 

0127 

0128 

0129 

0130 

0131 

0132 

0133 

0134 

0135 

0136 

0137 
0136 

0139 

0140 

0141 

0142 

0143 

0144 

0145 

0146 

0147 

0148 

0149 

0150 

0151 

0152 

0133 

0134 
0 ! j3 

0156 

0157 

0138 

0159 

0160 
0161 
0162 

0163 

0164 

0165 

0166 

0167 

0168 
0169 
C 1 70 
0 'I 71 

0172 

0173 

0174 

0175 

0176 

0177 

0178 


DO 30 J*1 ,NG 

30 WRITE (LUS , ' ( 13, F6 . 4 , 5 ( 1 PE1 1 . 3 ) > ' > 

1 J,ALP(J,1 ) ,U< J,t ) ,U(J,2> ,V< J,1 ) ,V<J,2> ,P< J) 

ENDIF 

ENDIF 

I F ( T . GT. THAX) GOTO 9999 

C TO SOLVE THE DIFFERENTIAL EQUATIONS 

C USING 2 STEP LAX-WENDROFF SCHEME 

C MacCornack's Method., BACKWARD PREDICTOR, FORWARD CORRECT 

C CENTER DIFFERFNCED ON TAU 

C FIRS: TO CET THE SPECIAL VARIABLES AND THEIR SPACIAL DERIVATIVES 

CALL DGCOEF < ALP , NG > ■ DRAG AND GENERIZED COFFF. 

CALL TAUQFW (MUEr , DR , R ,NG) ! STRESS 


VJT *VJ 

IF ( T .LT. TJ1 .OR, T . GT . TJ2) VJT*0 ! NO INJETCTIDN 
CALL JET ( BA , RHO , BR , V , VJT , NG) 

CALL FSOFW ( U , BA , R , NG ) 

C DETERMINE THE TIME-STEP SIZE 
IF ( NT .GT. 1) THEN 
CALL FN0DT ( DT DR , VDR , NG ) 

DT-DT+5 . #DTMIK/1 0 . *#NT 
IF ( DT .GE. D T MIN) THEN 
I*DLOG1 0 ( DT/DTM^N ) 

DDT*DTMIN«10.**I 
DT=*DINT ( DT/DDT+0 .001 ) *DDT 
IF ( DT .GT. DTMAX) DT-DTMAX 
ELSE 

WRITE(LUP , ' < 3X, A, 1 PEI 3 . 3 > ' > 

1 'STOP DUE TO TOO SMALL TIME STEP. DT*' , DT 

GOTO 9999 
ENDIF 
ENDIF 

C BACKWARD PREDICTOR 

CALL DWPDE ( DW , RDP , DR , DT , RH , NGM1 ) » INCREAMENT 

DO 40 J*1 , NGM1 
DO 40 1*1,5 

40 WP< J,I>»0 .5#<W< J + 1 , 1 )+U< J,I) )+DU( J,I ) * BASE-*- INCREAM 

C PREDICTION DATA COMPLETED , CONTINUE FOR CORRECTION 

CALL UOFW<WP,RH,MGMl ) 

CALL DGCOEF ( ALP , NGM1 ) 

CALL TAUOFU<MUEF,DR,RH,NGMl ) 

CALL JET ( BA , RHO, BRH, V, VJT, MG ) 

CALL FSOFW (UP , BA , RH , NGM1 ) 

C FORWARD CORRECTION 

CALL DWPDE < DW, P (2> , DR , DT, R (2) , NGM2) 

P<1)-PZERO 
DO 50 J-1,NGM2 


« MOMENTUM SOURCE 
» CONVECTIVE -F AND SOURCE-S 

! INITIAL INPLUSE TREAMENT 
• ROUND OFF TIME STEP 



i 
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50 


017? 

0180 

0181 

0182 

0183 

0184 C 

0185 


P ( J + l )*P< J) + < 0 . 23*<RDP ( J)+RDP ( J + l ) )+Q . 5*P < J>1 > > /* ( J ) 

DO 30 1-1,5 

WN< J + l , I . 23*< UP < J+l , I)+UP( J, I > ) + .5»<U< J + 1 , 1 > +DW ( J , I ) ) 
P(NG>-F<NCm i 

2ND STEP ( PREDICT I ON+C0RRECT ION ) COMPLETED 


018c 

0187 

3188 

018? 

0190 

01?1 

0192 
0173 
0 1 ?4 

0193 

0196 

0197 

0198 

0199 

0200 
0201 
0202 

0203 

0204 
02C3 
0206 

0207 

0208 

0209 

0210 
021 1 
0212 
021 3 
0214 
0213 
02!6 

0217 

0218 

0219 

0220 
0221 
0222 

0223 

0224 

0225 

0226 


IF< NT .Efl. 1 ) GOTO 10 

C Estimation of initial P completed, return to the initial condition 
C end start to advance The program in time. 

C DATA U AT THE NEW TIME STEP COMPLETED 

C IMPOSED B.C #6.4 

DU 1 * . 3* < UN ( 2 , 1 ) +W< 2 1 ) ) -0 . 

ALP <1 ,1 )*( . 5*( W< 1 , 3)+W(2,5) >-DWl *DT/DR >/RH< 1 ) 

IF ( ALP (1,1) .LT. ALMT(l)) ALP (1,1) -ALMT < 1 ) 

IF ( ALP (1,1) .GT. ALMT (2 ) ) ALP < 1 , 1) -ALMT ( 2 > 

UN < 1 , 5 ) -ALP ( 1 , 1 '*R < 1 ) 

DU 1-0 . -0 .5«(UN(NGM1 , 1 .'♦WCNGMl , | ) ) 

ALP ( NG , 1 )-< .5* (WCNGMl ,5>«-W(NG,5 ) >-DU1 »DT/DR )/RH(NGMl ) 

IF ( ALP ( NG , 1 ) .LT. ALMT < 1 > ) ALP ( NC , 1 ) -ALM f ( 1 ) 

IF< ALP (NG , i ) .GT. ALMT ( 2 ) ) ALP ( NC , 1 > -ALMT ( 2 ) 
k?NCNG,3)-ALP(NC, 1 )#R (NC) 

UN( 1 ,1 >-0. 1 NO RADIAL VEL . 

UNCI ,2)-0 . 

WN< NG , 1 ) -0 . 

UN<NC,2>-0 . 

WN( 1 , 3 )»UN< 2 , 3 ) *WN< 1 , 5 >/WN< 2 , 5 ) *R ( 1 )/R(2) 1 OMEGA-CON. 

UN ( 1 ,4)-UN<2, 4)*(R(1 >-UN(l ,3) ) ' ( R ( 2 ) -UN < 2 , 3 ) ) *R ( 1 >/W<2) 

IF < IU .EQ. 0) THEN 

UN < NG , 3 ) -WN( NGM1 ,3>»UN(NG,3)/WN<NGM1 , 5 ) *R ( NG > /R ( NGM1 ) 

UN(N£,4)i'Ut (NGM1 , 4 > » < R ( NG > -UN ( NC , 5) )/CR(NGM1 )-WN(NGHl ,3) ) 

1 *R(NG)/R (NGM1 ) 

ELSE 

UN ( NG , 3 > -9 . i NON-SLIP AT WALL 

UN(NG,4)»0. 

ENDIF 

C ARTIFICIAL TAMPING 
DO 60 1-1,5 

U(1 ,I)-<1 . -DAMP ) *WN( 1 , 1 )+DAMP*UN(2 , 1 ) 

U ( NG , I ) - ( 1 . -DAMP ) OWN < NG , I > ■►DAMP «UN ( NGM 1,1) 

DO 60 J»2,NGM1 

60 W< J, I )■< 1 . -DAMP )«UNt J , I )+DAMPo( WN( J-l , I ) +UN( J + l , I )-WN( J , I ) ) 


0227 

0228 

C 

SOLUTION FOR 

022? 

0230 

0231 

90 

T-T+DT 

0232 


GOTO 10 

0233 

0234 

999 

UK ITE ( LUP , 7 ) 

0233 


UR ITE ( LUP , 7 ) 

0236 


UR ITE ( LUP > 7 ) 

0237 

0238 

9999 

CONTINUE 


THIS TIME STEP COMPLETED 


• UPDATE TIME AND CONTINUE TO THE NEXT STEP 


' OPEN FILE FAILED ON FILE:' 
NAMR 

' IOSTAT- ' ,IOS 
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023 ? 

0240 

0241 

0242 

0243 

0244 
0243 

0246 

0247 
0246 
024? 
9230 

0231 

0232 

0253 

0254 

0233 
0256 
0237 
0256 
025? 
0260 
0261 
0262 

0263 

0264 

0265 

0266 

0267 

0268 
026? 

0270 

0271 

0272 

0273 

0274 
0273 

0276 

0277 
0270 
027? 
0260 
0261 
0262 

0263 

0264 
0283 
0286 
0287 
0266 
028? 
02?0 
0291 
02?2 

0293 

0294 

0295 

0296 

0297 

0298 


C CAUL £XEC( 1 1 » JTlHE , JTIME < 1 > ) 

C WRITE (LUS, ' (35X,3I4) ' ) JTIK€ 

CLOSE (LU3) 

END 


♦EMA /DAT A/ , /UWWW/ f /SOURCE/ 

SUBROUTINE INI! ,< 860423 . 1 537 > 

C TO SET-UP THE INITIAL CONDITIONS 

INTEGER I,IOS,ITLOC,IVTX,J,K,MH 
Z ,:TIH€,IW,LU?,LUS,NC,NCm 

PARAMETER (MM-101) 

CHARACTER NAMR*1£ 

REAL»8 AUNT, U,WP ,WN,DW ,RDP ,RH 
1 , U, V, ALP , P ,R , r,s 

4 , RUO , HUEF , V 1 8 , NA , ND , DO, GANNA 

7 , BR , BRH, DAMP, DR 

Y ,R J 1 ,R J2 , TJ1 ,TJ2,QJ,VJ 

X , A2 , 3ELA , DENI D2 ,D1 , D2 , P I , RE , RMIN , RTANK 

X ,D8,LS,V6,TS,PS 

X , ALP 1 0 , ONEGA , RPEAK , 9 JB , VPE AX 

X ,EVF( 2) ,VT,HU<2> ,HU1D2 


COMMON 

* /ALPLNT/ ALNT (2) 

Y /JETS/ RJl ,RJ2,TJ1 ,TJ2,QJ VJ 
Z /CONTP/ IU,LUP ,LllS,NG,NGHl # D'^NP,DR 

3 /DATA/ UtMM^) # V<HM,2>,ALP<MM,2> ,P(HH> ,*<MM> 

4 / DP AC 1 / RH0<4> ,HUEF(2> , Vl8<2> ,NA,ND 
3 /DSIZE/ DO ( 2 ) , GtWtMA ( 2 ) 

7 /SOURCE/ BR<HM) ,BRH(HH) 

9 /WWUW/ U(MH,S) ,WP(MM,3> ,WN(HH,5> ,DW<NN,3) , RDP < NH > ,RMiHM> 
DATA PI/3. 141596DB/ 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


>**•* NOTES: 
RHO 
ALNT 
DAMP 

DS ,LS , US ,TS 

RE 

RTANK 

IVTX 

RPEAK , VPEAK 
OMEGA 
EVF , HUEF 
DO , GAMMA 

RJl , RJ2 , TJl , 
TJ2,QJ,VJ 


PHASE- 1 »GAS , PHASE-2*LI0UID •**•** 

DENSITY, RHO<n < RHO< 2) 

LIMIT VALUES FOR ALP 1 , ALM ( 1 ) < ALP I < ALM< 2 ) 
NUMERICAL DAMPING FACTOR , NORMALLY -Q 
DEN91TY , LENGTH , VELOCITY AND TIME SCALES 
REYNOLDS #-V««RTANK»RMO< 2 > /HU < 2 > 

TANK RADIUS 

TYPE OF INITIAL FLOW. 

< 0»fi»»t , l*Pvr* r fit 4ti*n ,2-H . O . 3-GIT Vorft*) 

VORTEX PARAMETERS 

PURE ROTATION. V»OMCGA*r 

EFFECTIVE VISCOSITY , HUEF*< 1 ♦EVF ) «MU 

DIA. PARAMETERS: D-DO»ALP«»GAHMA 

TO DEFINE JET SIZE, PUMPING TIME, 

VOLUME FLOW RATE AND INJECTION MEAN SPEED 


5 FORMAT (2X,A,2<X,A) ) 

7 FORMAT (2X , A , 313) 

0 FORMAT (2X, A, 3(1 PE 12. 4) > 

9 FORMAT < A23 ,2<1PE13.4)> 

92 FORMAT ( X , 7 ( 1 PEI 1 .4) ) 
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~r 

* / 

i 


0249 

0300 C DEFINE THE PARAMETERS FOR THE PROBLEM. 

0301 

0302 RTANK-1.D0 • COULD BE SET TO 1 <m 

0303 RH0<2>-1 .0000+3 • * <KC/He#3> 

0304 MU<2)-1 .5140-3 ' * <KG/H-S> 

0305 

0306 C URITEd.8) 'ENTER RTANKOI) OR DEFAULT '.RTANK 

0307 C REAiH 1,1) RTANK 

0308 

030? C The »«lm of RTANK ,RH0<2> .end HU<2> ecu Id oil bo *•( to 1, oinco 

0310 C the length end density ocoloo oro booed on RTANK, end RHO<2> end the 

0311 C vole# of tho uisceslty MU<2) con bo conblnod into ond specified by tho 

0312 C Roynold* number RE. Thuu oil cherecter istic scsles(LS,VS,TS,e*«d DS) 

0313 C oro filed of tor RE is gluon. 

0314 

0313 RE-1.0D3 

0316 MIT£(1,8) 'on tor Roynol 

0317 REA0(1,e) RE 

0318 LS' "TANK 

0319 DS-RH0(2) 

0320 V3-RE«HU<2)/RH0(2)/LS 

0321 TS-LS/VS 

0322 P S«D8eVS**2 

0323 

0324 

0323 C A#TER THIS POINT AIL VARIABLES ARE BAESED ON THE CHARA. SCALES 

0326 C l.o. ALL VARIABLES ARE FINENSIONLE8S 

0327 

0328 

0329 DEN1D2-1 .29300/1 .000D3 • D1/D2 

0330 HU1D2-1 71D-3/1 ,314D- T • NU1/NU2 

0331 

0332 RHO(2)-RHO<2>/DS 

0333 M.‘(2)-NU(2)/<08nLSovS> • - 1 /RE 

0334 RHOd >-DENlD2eRMO<2> 

0333 HUd >-MUlD2eHU<2> 

0336 

0337 ALHTd >-0.00010 

0338 ALNT C 2 > *0 . 99999 

0339 D0d>-l.D-2/LS 

0340 DO<2)-1 .D-2/LS 

0341 GAMHAd >«2.D-1 

0342 GAMNA<2)-2.D-I 

0343 EVF (2>-l . 03 

0344 EVF ( 1 >-DENlD2/NUl D2*EVF<2> 

0343 DAHP-O.O 

0346 NA-4. 

0347 ND-4. 

0348 

9349 WRITE! 1 ,8) 'Enter DENSITY end VISCOSITY retie*.' 

0330 URIT£<1,8> 'D.F.-' ,9ENI02,MU1B2 

0331 READ! 1 ,e) 0CN1D2,NU1D2 

0332 

0333 URITEd.B) 'Enter BASE PiAMETERS: 001,002' 

0334 URITEO ,8) 'D.F.*', DO 

0333 REAOd.e) DO 

0336 

C337 URITEd.B) 'Enter SIZE EXPONENT: GAMMA! .6AHNA2' 

0338 Uei.Ed.B) '0. F. -'.GAMMA 


• MIN. OF ALP' 

• MAX. OF ALP 1 

• GAS DIAMETER et ALP l-l 

• LIQUID DIAMETER et ALP2-1 


• TURD. -PHASE-DISPERSION EFFECTS 
I MODEL 

• NUMERICAL DAMPING FACTOR (eg -.2) 

• UCICNTINC EXP. FOR ADM 

• HEIGHT INC EXP. FOR DRAG 


• LSee2/ ( NUeTS > -LSevS/NU 
e . , RE . D.F. -'.RE 

• LENGTH SCALE (M> 

• DENSITY SCALE (KC/M—3) 

• VELOCITY SCALE (H/S) 

• TIME SCALE <S> 

• PRESSURE SCALE 
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0359 READ<1,»> GAMMA 

0360 

0361 UR ITE (1,8) 'Enter weighting exponent: NA , HD , D . F . * ' , NA , ND 

0362 READ < 1 # * ) NA,ND 

0363 

0364 URITE(1,8) 'Ent er CAS VOLUME FRACTION 1 i *i T» ; ALMT 1 , ALMT2 . ' 

0365 URITE<1,8> 'D.F.*',ALMT 

0366 READ< 1 , • > ALMT 
9367 

0368 URITE<1,8> 'Enter eddy vimicoeity factor. D.F *',EVF 

0369 READ! 1 , * ) EVF 

0370 

0371 IU*1 

0372 URITE (1,7) 'Enter wall condi t ion , l*non*lip , 0**lip . D.F.',IU 

0373 README) IU 

0374 

0373 UR ITE< 1 .8) 'Enter nunorical danping factor. D.F.*', DAMP 

0376 READ* 1. *) DAMP 

0377 

0378 DO 10 K* 1 ,2 

0379 V18(K)-18.«MU(K) 

0380 10 MUEF ( K ) *HU < K > • < 1 . *E VK ( K ) ) * EFFECTIVE VISCOSITY FOR STRESS 

0381 

0382 RHO(3)-RHO( 1 )«RHO(2) 

0383 RM0(4)*RH0( 1 >-RHO<2) 

0384 

0383 RMlN-0.1 * MINIMUM FLOW RADIUS IN THE TANK 

0386 NG-101 ■ • OF GRID POINT* USED 

0387 NGM1*NC-1 

0388 DR*(1 . -RNIN)/NGM1 

0389 

0391 C Initial c leanning-op . 

C391 DO 13 J«1,MM 

0392 DO 15 K-l ,6 

0393 U( J , K )*0 . DO 

0394 13 U\ J ,K )«0 . DO 
0393 

0396 C MOMENTUM SOURCE, JET CONDITIONS 

0397 R Jl *8 . 5D- 1 

0398 RJ2-9.3D-1 

0399 VJ-10. • TANGENTIAL INJECTION SPEED 

0400 T J 1 *0 . 

0401 TJ2-10. 

0402 

9403 URI TE < 1 , 8 ) 'Enter JET SIZE defined by RJ1,RJ2. D.F.*' ,RJ1 ,RJ2 

0404 READ< 1 , * ) R Jl ,RJ2 

0403 URITE( 1 , 8 ) 'Enter INJECTION SPEED AND TIME RANGE, VJ,T1,T2' 

0406 U«ITE(1,B> DF.-' ,VJ,TJ1 ,TJ2 

0407 READd.o) VJ ; Tj1,TJ2 

0408 QJ»<RJ2-RJl )*UJ * JET VOLUME FLOW RATE 

0409 

0410 DO 20 J-t,NG 

0411 RC J>-RMIN*< J-l >*DR 

0412 RM<J>*R<J)*03«DR 

0413 BR< J)-DELA<<«< J> > ,DR,RJ1 ,RJ2>/<2.»PI> » JET DISTRIBUTION 

0414 20 BRH(J)*DELA<(RH<J>),DR,RJ1,RJ2)/(2.*PI) I PER RADIAN 

0413 

0416 

0417 C SETUP INITIAL CONDITIONS 

0418 IVTX«0 
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0419 

0420 

0421 

0422 

0423 

0424 
0423 

0426 

0427 

0428 

0429 

0430 

0431 

0432 

0433 

0434 
0433 

0436 

0437 

0438 

0439 

0440 

0441 

0442 

0443 

0444 
0443 
0446 
a447 

0448 

0449 
0433 

0431 

0432 

0433 
045' 
0453 
0456 
045 ^ 

0438 

0439 

0460 

0461 

0462 

0463 

0464 

0465 

0466 
P467 

0468 

0469 

0470 

0471 

0472 

0473 

0474 
0473 

0476 

0477 

0478 


OMEGA-0. 
VPEAK*Q . 
RPEAK-RNIN 


NANR-'Sinpln vortox' 

WRIT£CI,5) 'Ent«r data FILE NAHE for initial cond., l? any;' 
WRITE (1,5) ' D , F . ,NAHR 

READ ( 1 , ' < A ) ' ) NAMR 

IF ( NAMR .NE, .AND. NANR . N£ . 'Sinpl# vortor) THEN 

C INITIAL CONDITION FRON A GIVEN FILE NANR . 

OPEN ( 99 , FILE-NANR , I „3T IQS , STATUS* ' OLD ' , ERR=-299 > 

DO 23 J*1,NG I INITIAL VALUES FROM FILE NAUR 

23 StrtDC 99,-/ . . , ALP < J , 1 > > U ( J > 1 > > U ( J ,2) , V ( J , 1 ) ,V<J,2> 

CLOSE (99) 

ELSE 

C TO DEFINE INITIAL CONDITION. 

ALP1 0*2 . aO-l ‘ INITIAL GAS VOL. FRACTION 

WRITcM,8> 'Entor initial *lu# of alpl. D.F.*'.ALP10 
READ ( 1 , « ) ALP it< 

WRITE (1,7) 'Ent*r typo o7 vortox; 0*At r«*t,l-pur» rotation' 
WRITEC 1 , 7> '2-H.O. ,3-GXT. D . F. , IVTX 
READ( l,o) IVTX 
IF ( IVTX .GT. 0> THEN 
IF (IVTX .GT. 1> THEN 

WRITEC 1 ,8) 'En t«r PEAK TL and LOCATION for cla»mic vortoi' 
WRITE( 1 ,8) ' D . F . * ' , VPEAK , .. PEAK 
READCl,#) VPEAK ,RPEAK 

IF (RPEAK LE. 0.) RPEAK-1 . • SINGULAR AT ZERO 

ELSE 

WRITEC 1,8) 'En tor CIRCULAR SPEEDCr ad . /uni t r i *• ) . D . F . ■ ' 

1 , ONEGA 

READC 4,#) ONEGA 
END IF 
ENDIF 

DO 30 J»1 , NG 
ALP < J, 1 )*ALP 1 0 

VT-ONEGAoR 4 J) ■ PURE ROTATION 

IF C IVTX , GT, 0) THEN 
R JB-R ( J ) /RPEAK 

IF < IVTX .EG. t) THEN ! H.O. VORTEX 

VT-VT+l ,398*VPEAK/RJB*(1 ,-DEXPC-l . 23643*R JB**2) ) 

ELSE * G.I.T. VORTEX 

VT-VT+VPEAKOR JB»DEXP C ( 1 . -RJB**2)/2. ) 

ENDTF 

ENDIF 

DO 30 K*l,2 
IFCIW .EO. 0) THEN 

VCJ,K)*VT i NO WALL 

ELSE 

V(J,K)*VT*(1 ,-RCJ) )»*0. 1 ! BOUNDARY LAYER 

ENDIF 

30 U(J,K)-0. ! NO RADIAL VEL . 

ENDIF 

DO 40 J-1,NG * FORN W FOR NUNERICAL CAL. 
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r 


0 

i 

l 


•477 ALP(J,2>*1 .-ALP (J,l> 

•480 P(J)-0.D0 

•481 «(J,5)-*(J»«ALP(J,1) 

•482 DO 40 K-l ,2 

0483 *(J,K )-ALP( J,K)*U(J,K)*R< J) 

•484 40 W(J,K*2)-ALP(J,K>*V(J,K)«R( J> 

0483 

•486 

•487 C PRINTOUT PARAMETERS 

0488 WRITE (LUS, 3) 'INITIAL CONDITION PILE:', NAUR 

•487 WRITECLUS, 3> '••DIMENSION UNITS ARE IN MRS**' 

•490 WRIYE(LUS,7> 'DENSITY SCALE! kq/«3) ', DS 

•471 WRITE (LUS, 7) 'LENGTH SCALE -RTANR , («> ' ,LS 

•472 WRITE (LUS, 7) 'VELOCITY SCALE ( m /%) ' ,VS 

•473 WRITECLUS, 7) 'TIME SACLE(s>',TS 

•474 WRITE(LUS,7) 'PRESSURE SCALECPs) ', PS 

•473 WRITECLUS , •) 

•476 WRITECLUS, 7) 'R*yn«lds n»f*b*r, R*',RE 

•477 WRITE (LUS, 7) 'J*T nil*, RJ1 ,RJ2' ,RJ1 ,RJ2 

•478 MR I TEC LUS, 7) 'T*ft«*ftTi«i j*t, QJ, VJ' ,QJ,VJ 

•477 WRITE(LUS,7> 'Injsctlsn ti**,TJl,TJ2',TJl,TJ2 

•300 

•301 WR I TE C LUS , ' < /33X , "PHASE- 1 * , 8X , "PHASE-2* ) ' ) 

0302 WRITE (LUS, V ) ' Dtacx T y ' , RHOC 1 ) , RHOC2) 

•303 WRITE (LUS, 7) 'Viscosity ' ,mu 

0304 WRITE<LUS,7> 'Eddy visctsity f*ctor',EVF 

•3*5 WRITE(LUS,7) ' Ms* di« . ' , DO 

•316 WRITE (LUS, 7) 'Six* sip.', GAMMA 

•307 WRITECLUS.7) 'Ph«s* limits' , ALMT ( 1 ) , 1 . -ALMT (2) 

•308 WRITECLUS,*) 

•307 WRITECLUS, •> 

0310 WRITECLUS, '( " OTHER CONSTANS: IW, IVTX ,NA ,ND , DAMP , VPEAK , RPEAK " , 

•311 1 " , OMEC* , D1 /D2, MU1 /MU2" ) ' > 

•312 WRITECLUS, '(213,1 0,F3.2>'> IW, IVTX, NA,ND, DAMP 

•313 WRITECLUS, 72) VPEAR , RPEAK , OMEGA , DENI D2,MUl 02 

•314 WRITECLUS, ■> 

•313 UR I TE ( LUS , ' ( 1 1 • , F 1 • . 2 > ' > NG.RMIN 

•310 

•317 RETURN 

•318 

•317 277 WRITECLUP.S) 'OPEN FILE FAILED ON INPUT FILE:',NAMR 

•320 URITE(LUP,7> 'IOSTAT-' , IOS 

•321 STOP 1 1 1 

•322 END 

0323 

0324 *•*••**•*•••»*•••••••••*••••*•*«••••*•••***••••***•**•••*•*•**•***** 

•323 REAL *8 FUNCTION DELACR, DR, RJ1 ,RJ2> , (840423. 1337> 

•324 C TO DETERMINE THE EFFECTIVE NOZZLE SIZE AT EACH GRID LOCAT1 
•327 C THE SIZE IS IN THE FRACTION OF CRID SIZE DR (l.». KDELAC1 ) 

•328 

•327 REAL *8 R ,DR ,RJ1 , RJ2, R1,R2 

•330 

•331 R1-R-0.S*DR 

•332 R2-R1*DR 

•333 DELA-0 . D« 

4334 IFCRl CE UJ2 OR. R2 LE . RJ1 ) RETURN 

•333 IFCRl . L T . RJ1 ) Rl-RJl 

•334 IFCR2 CT. RJ2) R2-RJ2 

•337 DELA-(R2-R1 )/DR 

•338 RETURN 
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0539 

0540 

0541 

0542 

0543 

0544 

0545 

0546 

0547 

0548 

0549 

0550 

0551 

0552 

0553 

0554 

0555 

0556 

0557 

0558 

0560 

0561 

0562 

0563 

0564 

0565 

0566 

0567 

0568 

0569 

0570 

0571 

0572 

0573 

0574 

0575 

0576 

0577 
0570 

0579 

0580 

0581 

0582 

0583 

0584 

0585 

0586 

0587 

0588 

0589 

0590 

0591 
0392 

0593 

0594 
0593 
0596 
0397 
0598 


2ND 


*******««*«***#»###**#* «**»««»«*#**#**?*»««#* »#**»**«***«*#**** 
SUBROUTINE DERIV1 <Y,DY,DX,N2) , <860425, 1337> 

C GET 1ST DERIVATIVE, USING CENTERED DIFFERENCE 
REAL#0 DX,Y<1),DY(1),C 
ENA Y , DY 

C*5 . D- 1 /DX 
DO 10 J*2,N2-1 

10 DYC J)«C*< Y(JM)-Y(M ) ) 

*V CALL DUSUB< Y <3),1,Y,1 , DY(2) ,1 >N2-2) 

DY<1)*(Y(2)-Y(1))/DX ' BASED CN 3-END PTS 

DY(N2)*(Y( N2) -Y f N2-1 ) )/DX 
DYU >-2.#DY<l >-DY<2) 

DY(N2)»2. #DY<N2)-DY<N2-1 ) 

*V CALL DWSHY<5.D~1/DX,DY,1 ,DY,1 ,N2) 

RETURN 

END 


•ENA /COEFF/ 

SUBROUTINE DGCOEF < ALP , N2) , <860423. 1537> 

C CALCULATE THE DRAG, ADDED HASS AND SENEJIZED COEFF, 

INTEGER J , MM , N2 
PARAMETER <MM»101 ) 

REAL#8 ALP ( MM , 2 ) 

EMA ALP 

REAL*8 C, CP A, CD ; RHO , MUEF , VI 8 , NA , ND 
X ; AA ; AA1 ,AA2,AD,AD1 , AD2 , A1 , A2 , At 2 , DAI 2 ; DB2 , D1 ,D2,WT1 ,WT2,X 

COMMON 

1 /COEFF/ C(MM,2,2) ;CPA(MH,2) ;CD(MM;2) 

4 /DRAG1/ RH0C4) ;MUEF<2) , V18<2) ,NA,ND 

DO 50 J»1;L2 
A1»ALP(J; 1 ) 

A2-ALP( J;2) 

A12-A1«A2 

C TO GE T DRAG COEFF. AD 
CALL SIZES<D1 ,D2 ; A1 > 

AD-V18<2)»A3/<A2#D1*D1 ) f -ADI IF A2>.78 

IF( A3 .LT. .78) THEN 
AD2-V18U >«A2/< (1 . -A2/ , 8) *D2 > »#2 
2F( AD2 .LT. AD) THEN 
UT1-A2««ND 
UT2-A1*»ND 

AD- < AD»UT 1 ♦ AD2«UT2 > / < WT 1 +W T2 > 

ENDIF 

FNDIF 

C ADDED MASS COEFF. AA 
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If 


0599 

0600 
0601 
0602 

0603 

0604 

0605 

0606 

0607 

0608 

0609 

0610 
0611 
0612 

0613 

0614 

0615 

0616 

0617 

0618 

0619 

0620 
0621 
0622 

0623 

0624 

0625 

0626 

0627 

0628 

0629 

0630 
o*3: 

0632 

0633 

0634 

0635 

0636 

0637 

0638 

0639 

0640 

0641 

0642 

0643 

0644 
06*5 

0646 

0647 

0648 

0649 

0650 

0651 

0652 

0653 

0654 

0655 

0656 

0657 

0658 


AAl*Al2»RH0(2>/<Al*A2/< • 5«*1 5*A1 ) ) 
AA2»A12»RH0<1 )/CAl/( .541 .5*A2)*A2) 

WT1*A2**MA 

UT2»A1««NA 

AA« ( AA 1 »WT 1 +AA2»W T2 > / < WT 1 +UT2 ) 

C THE GENERIZED COEFF. CP A, C, AND CD 

DB2*A12«RH0<3>*AA*<RH0( 1 )»A1*RH0(2>«A2> 
CP A < J , 1 )-Al*<A12*RH0(2)>AA)/DB2 
CPA< J,2>-A2»<A12*RH0<1 )*AA)/DB2 
C<J,1 ,2>*AA/DB2 
C<J,2,1)*C(J,1,?> 

C<J,1 ,1 >«A2*RH0<2)/DB2*C(J,l i 2> 

C<J ,2,2)»Ai**H0n)/VB2+C<J , 2,1 > 

CD < J , 1 > »-A2*RH0 ( 2 ) *AD/DB2 
CD < J , 2 ) -A 1 #RHO < 1 ) « AO/ DB2 

50 CONTINUE 

RETURN 

END 


♦EHA /DATA/ , /FANDS/ , /TAU/ , /COEFF/ 

SUBROUTINE DWPDE ( DU , RDP , D£ , DT , RR , N2) ,< 860425 . 1 537 > 

C 

C TO GET DW OF THE PDE* 

C 

INTEGER J, JP1 ,K,KP2,MM,N2 
PARAMETER <MM»101> 

k ^ ' *8 DW<MM,5> ,RDP(NM> , PR,DT,RR<MM> 

ENA DW > RDP / RR 

C NOTES: COEFF. C *C»ALP WHEN THIS IS CALLED 

REAL*8 C, CPA, CD, U,V,ALP,P,R, F,S 
8 , TRR , TRA, TAA,RTRR , RTRA 

4 , RHO , NUEF , V 1 8 , NA , ND 

X ,ALP1 ,ALP2,CP1 ,CP2,0TDR,DU1 ,DW1DT,G1 , C2 ,MDT ,UT,MJ1 , UJ3,WJ4 


COMMON 

1 /COEFF/ C<MM,2,2),CPA<HH,2>,CD<MM,2> 

3 /DATA/ U<HH,2),V<MN,2>,ALP<MM,2>,P<HH>,R<HM) 

4 /DRAG1 / RH0<4),HUEF<2> VM8<2),NA,ND 
6 /FANDS/ F ( MM , 5 ) , S( MM , 5 ) 

8 /TAU/ 7kR < MM , 2 ) , T R A < MM , 2 > , TAA< MM , 2 ) , RTRR < MM , 2 > , RTR A < MM # 2 ) 

DTDR-DT/DR 

HDT-0.5*Dr 

DO 10 J»t,N2+1 • CHANGE C TO ALP»C 

C< J, 1 , 1 )*ALP < J, 1 )*C( J, 1 , 1 ) 

C ( J , 2 , 1 ) "ALP <J,2)4C(J,2>1 > 

C< J, 1 ,2)-ALP< J, 1 >*C<J,* ,2) 

1 0 C < J , 2 , 2 ) "ALP < J , 2 > *C < J ,2,2) 

DO 20 J-1 , N2 
JPW+1 

DW ( J , 5 ) »DTDR* < **F ( JP 1 , 5 > *F ( J , 5 > ) +MDT# < S < JP 1 ,5)*S<J,3> ) 

DO 25 K*1 , 2 
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0659 

0660 
0661 
0662 

0663 

0664 

0665 
Q666 

0667 

0668 

0669 

0670 

0671 

0672 

0673 

0674 

0675 

0676 

0677 

0678 

0679 

0680 
0681 
0662 

0683 

0684 

0685 

0686 

0687 

0688 

0689 

0690 

0691 

0692 

0693 

0694 

0695 

0696 

0697 

0698 

0699 

0700 

0701 

0702 

0703 

0704 

0705 

0706 

0707 

0708 

0709 

0710 

0711 

0712 

0713 

0714 

0715 

0716 

0717 
0710 


G1-.5*(C< JP1 ,K,1 )+C(J,K,1 )> 

G2».5*<C< JP1 >K t 2)+C (J,K,2>) 

DU< J,K>*DTDR#(-F< JP1 , K >+F ( J , K ) +G1 »< RT RR ( JP 1 , 1 ) -R TRR ( J , 1 ) ) 

1 +G2* ( RTRR ( JP 1 ,2>-RTRR( J,2) )) 

2 +HDT»<S(JP1 ,K)+SCJ,K>) 

KP2-K+2 

25 PW< J,KP2)»DTDR*<-F< JP1 , KP2 > +F ( J , K P2 ) 

1 +G1*(RTRA< JP1 , 1 >-RTRA< J, 1 > > 

2 ♦G2«(RTRA< JP1 , 2) -RTRA< J , 2 ) ) ) 

3 ♦HOT*(S< JTP1 ,KP2>+S< J,KP2) 7 

C DP FOR PRESSURE CORRECTION 

CP1-0.5*(CPA< J, 1 ) +CP A< JP1 , 1 ) ) 

CP2-0.5*(CPA<J,2>+CPA<jri ,2) ) 

IF < -DU<J,:> .GT. DU<J,2)> DU( J, 1 )»-DW< J,2 ) ! DP>*0 

RDP< J)«<DW< J, 1 )+DW< J,2> )/CCPl+CP2) 

DW ( J , 1 ) *DW ( J , 1 ) -CP 1 *RDP < J ) 

DW ( J , 2 ) a -DW< J , 1 ) 

20 RDP< J>-RDP(J)/DT 

RETURN 

END 


*##«***#*#*##*«*#•*#**#*#******** »«*#*****###*###****#*********#*##* 
♦ENA /COEFF/,/pATA/./WWWW/ 

SUBROUTINE FNDDT < DT , DR , NDR , NG ) , < 860 425 . 1 537 > 

C DETERMINE THE TIME-STEP SIZE 


INTEGER I, J,LUP,MM,NC 
PARAMETER (MM* 101) 

REAL»8 DT,DR,NDR<2> 

REALMS C, CPA, CD, RHO , MUEF , V 1 8 , NA , ND 
9 , U , UP , UN , DU, RDP , RH , U,V,ALP,P,R 

X , DUM1 , DUM2 


COMMON 

2 /COE FF/ C<MM,2,2> ,CPA<MM,2) ,CD<MM,2> 

3 /DATA/ U(MM,2) / U(MM,2) / ALP<MM,2),P<MM) ,R<MM) 

4 /DRAGl / RHO< 4) ,MUEF<2) , VI 8 ( 2 ) , NA ,ND 

9 /UUUU/ U ( MM , 5 > ,UP ( MM , 5 ) , UN< MM t 5 ) , DU< MM , 5 ) , RDP < MM ) , RH( MM ) 
DATA LUP/1/ 


DUM1-0. 

DO 10 J-1,NG 
DUM2— CD( J,1 ) 

1 +DABS ( U ( J , 1 ) ) /DR 

2 ♦UDR < 1 # * ALP <J,1)*C(J,1,1) 

3 ♦NDR ( 2 ) »ALP ( J , 2 ) *C < J , 1 ,2) 
IF < DUM1 .LT. DUM2) DUM1-DUM2 
DUM2-CD ( J , 2 > +DABS < U < J , 2 ) > / DP 

1 ♦MDR ( 2 > *ALP <J,2)*C(J,2,2) 

2 ♦000(1 )*ALP< J, 1 >»C< J,2, 1 ) 
IF ( DUM1 . LT . DUM2 ) DUM1-DUM2 

10 CONTINUE 


< t/St 

i CONNECTIVE 
' NDR»2*MUEF/DR**2 
• CROSS VISICOSITY 


C FIND THE MAXIMUN OF DU 
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VECTOR OPERATION 


0719 

0720 

0721 

0722 

0723 

0724 
0723 

0726 

0727 

0728 

0729 

0730 

0731 

0732 

0733 

0734 

0735 

0736 

0737 

0738 

0739 

0740 

0741 

0742 

0743 

0744 
0743 

0746 

0747 

0748 

0749 

0730 

0731 

0732 

0733 

0734 
0753 

0736 

0737 
0758 
0739 

0760 

0761 

0762 

0763 

0764 

0765 

0766 

0767 
0760 

0769 

0770 

0771 

0772 

0773 

0774 
0773 

0776 

0777 

0778 


CALL DWMAX<I,DW,1 ,NC> 

»V DUM 1 *D ABS ( DU < I , 1 ) ) 

*V CALL DWMAX ( J , DU ( 1 , 2 ) , 1 , NG ) 
* V DUM2*D ABS ( 04 < J , 2 ) ) 

DT-1 . D0/DUM1 

RETURN 

END 


♦ENA /DATA/ , /FANDS/ , /TAU/ , /CP C FF/ 

SUBROUTINE FSOFW< W , BA , Rk , N2 > ,< 060423 . 1 537 > 

C 

C CALCULATE THE CONVECTIVE-F AND SOURCE-S TERMS 

C 

INTEGER J,K,KP2,L,MM,N2 
PARAMETER (MM-101) 

REALMS U,V,ALP,P,R, F,S, RHO , MUEF , VI 8 ,WA ,ND 

2 ,TRR,TRA,TAA,RTRR,RTRA, C, CPA. CD 

X ,ALPD<2>,RDU,RDV 

REALMS W < MM , 5 ) , BA < MM , 2 ) , R R < 1 > 

EMA W,BA,RR 

COMMON 

1 /COEFF/ CCMH,2,2> ,CPA(MM,2) ,CD<MM,2) 

3 /DATA/ UCMM^ljVCMM^^ALPCMM^^PCMM^RCMM) 

4 /DRAG1/ RH0<4) ,MUEF<2> , V18C2) ,NA,ND 
6 /FANDS/ F ( MM , 5 ) , S < MM > 3 ) 

8 /TAU/ TRR (MM,2> , TR A < MM , 2 ) , TAA < MM , 2 ) , R TRR ( MM , 2 ) , RTR A < MM , 2 ) 

*V CALL DWMPY(W,l ,U,1 ,F , 1 ,2*MM) 

*V CALL DWMPY < W( 1 ; 3) , 1 , U , 1 ,F< 1 , 3) , 1 , 2*MM> 

*V CALL DWMO V <U,1,F<1,5),1,MM) 

DO 20 J«1 ,N2 
F( J,3)*W< J, 1 ) 

S(J,5)*0. 

RDU«RR<J)#<U< J, 1 >-U< J,2> ) 

RDV-RR(J)»<V(J, 1 >-V< J,2) > 

DO 20 K*1 ,2 
KP2*K+2 

F< J,K)-W<J,K)#U< J,K) 

F( J,KP2)«U< J,KP2>»U< J,K > 

S(J,K)»ALP ( J,K >*<V< J,K )##2+CD< J,K >*RDU-C< J,K , 1 >»fAA< J, 1 > 

1 -C( J,.<,2)*TAA<J,2>> 

20 S(J,KP2)-ALP< J,K )#<-U< J , K )*V( J ,K )+CO< J , K > #RDV+C ( J , K >*<BA<J,1 ) 

1 +TRA< J, 1 ) >+C<J,K,2>»(BA(J,2>+TRA<J,2>) > 


RETURN 

END 


SUBROUTINE JET < BA # RHO , BR , V , VJ , NG ) 
C 

C INJECTION MOMENTUM SOURCE 

C 
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0779 

0780 
0701 
0782 

0703 

0704 

0705 

0706 

0787 

0788 

0789 

0790 

0791 

0792 

0793 

0794 
0793 

0796 

0797 

0798 

0799 
0600 
0801 
0802 

0803 

0804 

0805 

0806 

0807 

0808 

0809 

0810 
0811 
0812 

0813 

0814 
0015 
0316 
0817 
0816 

0819 

0820 
0821 
0822 

0823 

0824 
0823 
0826 

0827 

0828 
0029 

0830 

0831 

0832 

0833 
0034 
0833 
0836 
0037 
0838 


PARAMETER <MM*1B1) 

REAL* 8 BA < MM, 2) , RHO< 1 ) , Bfc \ 1 ) , V ( MM , 2 > , V J 

EMA BA,BR,V 

do to j-i 

I F< BR ( J ) .GT, 0.) THEN 

Q«BR<J>*VJ ' Plow r ort/radion 

BA< J, 1 )*0«< VJ-V< J, 1 ) )*RHOC 1 ) 1 NET MOMENTUM GAIN 

BA<J,2)-Q»CVJ-V< J,2> )*RH0(2) 

ELSE 

BA< J,1 )*0.D0 
BA<J,2)-0.DO 
END IF 

10 CONTINUE 
RETURN 
END 


****»#*##********##*#****#**#c*#« *************** 

SUBROUTINE SIZES < D1 , D2 , ALP 1 ) , <860423. 1337) 
C TO DETERMINE THE PARTICLE DIAMETERS 

C 

REAL*8 D1 > D2, ALPt , DO , GAMMA 
COMMON /DSIZE/ DO (2) ,GAMMA<2) 

Dl*DO< 1 ) *ALP 1 **GAMMA< 1 ) 

D2*D0 ( 2 ) * < 1 . 0D0-ALP1 >**GAMMA<2> 

RETURN 

END 


♦EMA /DATA/ , /TAU/ 

SUBROUTINE TAUOFW< MU , DR , RR > N2 ) , < 860425 . 1 337 > 
C STRESSES AND THEIR DERIVATIVES 

C 

REAL*8 MU(2) ,DR,RR<1 ) 

EMA RR 

PARAMETER CHH-101) 

REAL*0 U,V,ALP,P,R 
2 ,TRR,TRA,TAA,RTRR,RTRA 

X , TAMU, TMU 


COMMON 

3 /DATA/ V < MM , 2 ) , V < MM , 2 ) . ALP < MM , 2 ) > P < MM ) , R ( MM ) 

8 /TAU/ TRR < MM, 2) ,TRA<HM,2> ,TAA<Mh,2) ,RTRR (MM, 2) ,RTRA<MM,2) 


DO 50 K»1 ,2 

CALL DERI VI ( U < 1 , K ) , TRR < 1 , K ) , DR , N2> 

*V CALL DUMP Y < ALP ( 1 ,K ) , 1 * TRR < 1 , K > , 1 , TRR ( 1 , K ) , l , N2 ) 

*V CALL DWSHY<2 . *MU<K ) ,TRR< 1 ,K ) , 1 ,TRR< 1 ,K ) , 1 ,N2) 

•V CALL DWDIV<V<l,K>,l,RR,l,rAA<l,K>,i,N2>' V/R 

DO 10 J»1 ,N2 

1 0 T A A ( J , K ) ■ V ( J , K > /RR < J > 

CALL DERIV1 (TAA< i ,K > , TRAC 1 ,K>,DR,N2> 


' dU/dR 
* *ALP 

i »2*MU*TRR 


i d<>/dR 


TMU*2 . *MU<X ) 

DO 20 J*1 , N2 
TAMU»TMU*ALP( J,K) 

TRR < J , K ) «TAMU»TRR < J , K ) 
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i 


t 


083? 

0840 

0841 

0842 

0843 

0844 

0845 

0846 

0847 

0848 

0849 

0850 

0851 
0352 
065? 
0854 
0853 

0856 

0857 

0858 
0839 
0860 
0861 
0862 

0863 

0864 
0863 
0866 

0867 

0868 

0869 

0870 

0871 

0872 

0873 

0874 

0875 

0876 

0877 
0678 
0079 
0680 
0881 
0882 

0883 

0884 

0885 

0886 
0887 
0868 

0889 

0890 
0091 

0892 

0893 

0894 
0893 
0876 
0897 
0398 


TRAC J,K >-MUCK ?#ALPC J,K )*RRC J)»TRAC J,K ) 

T AA C J , K ) « TAMU«U C J , K > /R R < J ) 

RTRRC J,K ) »RR (J)#TRR(J,K) 

20 RTRAC J,K )-RRCJ>*TRAC J,K> 

i *R 
• *ALP 
' *MU*TR A 


f #ALP 
» »2*MU-TAA 

i R*TRR 
» R*TK A 

50 CONTINUE 
RETURN 
END 


*V CALL DWMPY<RR,1, TRAM, K ) , i , TRA < 1 , K ) , 1 , N£) 

#V CALL DUMPY ( ALP ( 1 ,K> , 1 , TRAC 1 # K) ,1 , TRACI ,X) , 1 ,N2> 

*V CALL DUSMYCMUOO , TRAC 1 ,K > , 1 ,TRAC 1 ,K ) , 1 ,N2) 

*V CALL DUDIVCUC 1 ,K ) , 1 ,RR , 1 , TAA( 1 ,K ) , 1 ,N2) ' U/P 

•V CALL DUMP Y C ALP < 1 , K ) , 1 , TAAC 1 , K ) , 1 , TAAC 1 , X ) , 1 , N2) 

»V CALL DUSMYf 2. #HU(K ) ,TAA< 1 ,K ) , 1 ,TAAC 1 ,K ) , t ,N2> 

»V CALL DUMPY C RR , l , TRR C 1 , K ) , 1 , RTRR < 1 , K ) , 1 , N2 ) 

*V CALL DUMPYCRR, 1 ,TRAC 1 ,K ) , 1 ,RTRA< 1 ,K ) , I ,N2> 


♦EMA /DATA/ 

SUBROUTINE UOPUC M, RR > N2) , C 860425 . 1 537> 

C CONVERTS U TO THE INDEPENDENT VARIABLES ( U , V , ALP ) 

PARAMETER CMM»101) 

REAL»8 UCMM,5> ,RRCMM) 

EMA U ,RR 

REAL*8 ALMT, U,V,ALP,P,R 
X t U J 6 

COMMON 

X /ALPLMT/ ALMTC2) 

3 /DATA/ U ( MM , 2 ) , V C MM , 2 > , ALP C MM > 2 ) , P < MM ) , R C MM ) 

•V CALL DUDIVCUC 1 ,5) , 1 ,RR , 1 , ALP # 1 , N2 ) 

C CHECK VOLUME FRACTION 6 FLOU DIRECTIONS 

DO 30 J-I , N2 
ALPC J, 1 )-Uc\: j)/RR C J ) 

IF C ALP C J , 1 J .LE. ALMTC1) .OR. ALPCJ,1> .Cl. ALMT C 2 ) ) THEN 
IF C ALP C J , 1 ) .LE. ALMTCl)) THEN 
ALP C J , 1 >«ALMT< 1 ) 

UC J,3)»UCJ,4)*ALMTC 1 )/C 1 . D0-ALMT C 1 ) ) 

ELSE 

ALPCJ, 1 J-ALMTC2) 

UC J.AJ-UC J,3)*C 1 . DO-ALMT ( 2 ) )/ALMTC2) 

END IF 

UC J,1 )>0.00 
UC J , 2 ) ®0 . DO 
UC J,3>»ALP< J,1 )*RRC J) 

ENDIF 

ALP < J ,2 ) * t . DO-ALP < J , 1 ) 

IF C WU,2> LT. 0.) THEN • PHASE-2 DOES NOT MOVE IN 

UC J t 2)-0. 

UC J > t >«0 . 

ENDI r 
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r 


089? 


U< J, 1 )»U< J, 1 )/W( J,5) 

0900 


v<jr, 1 )-U< J , 3 > /W< J , 3 > 

0901 


UJ6»RR< J)#ALP< J,2> 

0902 


U(J,2)*VlJ,2)/UJt> 

0903 


J,2)-W< J,4)/W.T6 

0904 

50 

CONTINUE 

0903 



0906 


CAuL DWDIWW, ; ,y<l ,5) , 1 ,U,1 ,N2) 

0907 

*v 

CALL DUDIV<Wn , 3) , 1 , W< 1 , J ) , 1 , V, 1 , N2) 

0908 


CALL DUMP Y ( ALP < 1 ,2),1,RR,1,W<1,6),1,N2) 

0909 


CALL 0W0IW(a(l,2) .1 ,U<l,6>,t,U<l,2>,l,N2> 

0910 


CALL DUDIV(W<1 ,4>,t ,U(1 ,<>>,1 ,UC1 ,2>,1 ,N2> 

091 1 



0912 


RETURN 

0913 


END 

0914 
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Exhibit A 


A Sample Input 


:GLVM 

Enter la 'or winq data. D.F.= b 
SO 

Enter r ILE NAME tor saving aa ta. 

TS 153:: LB 

Enter NGTES<<73 CHAR.) for the job 
SAMPLE RUN OF TEST 153- 

enter Reynolds no., RE. D.F.= 1.0000L+05 

Enter DENSITY ana VISCOSITY ratios. 

D.F.- 1.2930E-03 1.1295E-02 

’’Enter BASE DIAMETERS: D01.D02 
D.F.- 1.0000E-02 1.0000E-02 

’’Enter SIZE EXPONENT: GAMMA 1 , GAMMH2 
D. F. = 2 . 00C0E-01 2.0000E-P1 

• * 

Enter weighting exponent: NA.ND. D.F.= 4.0000E+00 4 . 0000E+00 

’’Enter GAS VOLUME FRACTION limits : ALMT 1 . ALMT2 . 

D.F.- 1.0000E-04 9.9999E-01 

» * 

Enter eddy visicosity factor. D.F.* 1.!443t+'.i2 1.0000E+03 

1000.1000 

Enter wail condi tion ,.l -nonsl ip . 0 = si ip . D.r. 1 

• 

Enter numerical damping factor. D.F.* O.OQOOt+OO 

Enter JET SIZE defined by RJ1.RJ2. D.F.- 3.5000E-01 3.5000E-01 

’’enter INJECTION SPEED AND TIME RANGE. Vj.Tl.T2 
D.F.- 1 . OOOOE+GI O.OOOOE+OO 1.0000E+01 

1 . 0.1 

Enter data FILE NAME for initial cond . i ; any: 

D.F.= Simple vortex 

Enter initial vaiue of aipl . D.F.* 2.5000E-01 

» 

Ente r type of vortex: 0=At rest. Impure rotation 
2-H.O. ,3-GIT. D.F.- 0 

Enter INITIAL and FINAL TIMEs. D.F.* O.OOOOE+OO 5.00U0E+00 

0..01 

Enter TIME STEP for output. D.F.= 2.0000E-01 

.01 

I NP= 1 T-0.000E+00 NT- 2 DTM . C00E-06 

1 NP- 2 T»1 .200E-02 NT- 6 DT-3.000E-03 
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Exhibit B 


f 

* 

i t 



A Sample Output 


TS 153 

T -00004 IC UN CR LB US INC 

00024 RLKS R = ' 

j o r, •'! 

0001 

1 SAMPLE RUN OF TEST 153 



0002 

INITIAL CONDITION FIlE: . 



0003 

••DIMENSION UNITS ARE IN MKS*« 


0004 

DENSITY SCALE(kg/m3) 

1 .OOflOE + 03 


0005 

LENGTH SCALE-RTANK.(m) 

1 . 0000E+00 


0008 

VELOCITY SCALE < m/s) 

1 . 5 1 40E-0 1 


0007 

TIME SACLE(s) 

6 .6050E <-00 


0008 

PRESSURE SCALE < Pa > 

2.2922E+C 1 


0009 




0010 

Reynoids number. Re 

? .OOOOE + 05 


0011 

Jet size, RJ1 ,RJ2 

8.5000E-01 

3.5000E-01 

0012 

Tangential jet, QJ.VJ 

1 . OOOOE-0 1 

1 .0000E+00 

0013 

Injection time,TJ1,TJ2 

0 . OOOOE+OO 

i .noooE+oo 

00i4 




0015 


PHASE- 1 

PHASE-2 

0016 

Densi ty 

1 . 2930E-03 

1 . 0000E+00 

0017 

V iscos 1 ty 

1 . I295E-07 

! . 9000E-05 

0018 

Eddy viscosity factor 

1 . Q00QE+03 

1 . 0000E+03 

0019 

Base dia. 

1 . 0000E-02 

! . 0000E- 02 

002C 

Size exp. 

2.0000E-01 

2 . 0000E-01 

0021 

Phase limits 

! . 0000E - 04 

1 . 00 1 4E-05 


t'Ucc 

0023 

0024 OTHER CON? TANS: IW. I VTX ,NA,ND, DAMP . VPEAK .RPE4K , OMEGA . D1 /D2 ,MU ! /MU 


0025 

0026 

1 0 4. 

0 . OOOOE+OO 

4. 0.00 
I . OOOOE-O 1 

0. OOOOE+OO 

I . 2930E-03 

1 .1295E-02 


0027 

0028 
0029 

1 NP* 

101 

.10 

T *0 . 000E+00 

NT* 2 

DT*1 .000E- 

06 


0030 

J 

ALP1 

U1 

U2 

Vi 

V2 

p 

0031 

I 

.2500 

0 . 000E+00 

0 .000E *00 

0.000E+00 

0 . OOOE+OO 

0. 000E+00 

0032 

2 

.2500 

0 - 000E+00 

0.000E+00 

0.000E+00 

0.000E+00 

0. OOOE+OO 

0230 

97 

.2500 

OUTPIJ 

-8.685E-10 

U IN THE BETWEEN OMITTED 

2.895E-10 3.513E-04 3.474E-04 

3.969E-05 

0231 

98 

. 2500 

-t . 775E- 1 0 

5.918E-1 1 

1 . 268E-04 

l . 250E-04 

3.972E-05 

0232 

99 

.2500 

-2.785E-11 

9.284E-12 

4.319E-05 

4.256E-05 

3.973E-05 

0233 

100 

.2500 

-5.225E- 12 

1 .742E-12 

1 .553E-05 

1 .53 IE-05 

5.973E-05 

C234 

101 

.2500 

0 . 000E+00 

0.000E+00 

0 . 000E+00 

0. 000E+00 

3.973E-A5 


THE REST OF Th'E OUTPUT IS OMITTED 
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Complete step 
Correction step 

Prediction step 


Initial value 


Flgur« 1. TWo Stop Dlfforcnco Schoet (Backward Prodictor - 
Forward Corractor Varalon) 
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I 
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Tli« annular raglon batwean two daahad drdas is 
tha raglon of injactlon. 



‘UOIjODJJ 9UJn|0A SDQ 
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Figure 3. Gas Voluae Fraction Distributions 



